Collective modes and superflow instabilities of strongly correlated Fermi superfluids 
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We study the superfluid phase of the one-band attractive Hubbard model of fermions as a pro- 
totype of a strongly correlated s-wave fermion superfluid on a lattice. We show that the collective 
mode spectrum of this superfluid exhibits, in addition to the long wavelength sound mode, a sharp 
roton mode over a wide range of densities and interaction strengths. We compute the sound velocity 
and the roton gap within a generalized random phase approximation (GRPA) and show that the 
GRPA results are in good agreement, at strong coupling, with a spin wave analysis of the appro- 
priate strong coupling pseudospin model. We also investigate, using this two-pronged approach, 
the breakdown of superfluidity in the presence of a supercurrent. We find that the superflow can 
break down at a critical flow momentum via several distinct mechanisms — depairing. Landau in- 
stabilities or dynamical instabilities — depending on the dimension, the interaction strength and 
the fermion density. The most interesting of these instabilities is a charge modulation dynamical 
instability which is distinct from previously studied dynamical instabilities of Bose superfluids. The 
charge order associated with this instability can be of two types: (i) a commensurate checkerboard 
modulation driven by softening of the roton mode at the Brillouin zone corner, or, (ii) an incommen- 
surate density modulation arising from superflow-induced finite momentum pairing of Bogoliubov 
quasiparticles. We elucidate the dynamical phase diagram showing the critical flow momentum of 
the leading instability over a wide range of fermion densities and interaction strengths and point out 
implications of our results for experiments on cold atom fermion superfluids in an optical lattice. 



I. INTRODUCTION 

The study of strongly correlated fermionic superfluids 
is of interest in the context of solid state materials, such 
as cupratei and pnictide superconductors,- as well as ul- 
tracold atomic gases in optical lattices.'^''* One of the key 
goals of condensed matter physics is to understand the 
nature of single particle and collective excitations in such 
strongly correlated systems. Another important problem 
is to eludicate the mechanisms by which superfluidity 
breaks down in these lattice systems — such an under- 
standing would shed light on the critical current and on 
the nature of vortex cores in strongly interacting super- 
fluids,- which are issues of significant experimental and 
theoretical interest. Moreover, the inverse critical su- 
perflow momentum can be shown to be directly related 
to characteristic length scales, associated with the low- 
energy excitations in the system and can thus serve as 
a useful probe of these excitations. This is especially 
relevant for neutral cold atom superfluids, in which the 
critical flow is unaffected by complications such as disor- 
der effects and current induced magnetic fields which are 
present in solid state superconductors. While there has 
been a significant amount of experimental and theoretical 
progress in understanding the breakdown of superfiow for 
Bose superfluids in an optical lattice r'liiSii^ and for Fermi 
superfluids in the absence of a lattice potential )^"'^^ this 
issue has not been addressed in detail in the context of 
Fermi superfluids on a lattice. 

In this paper we study the collective modes and the 
mechanisms for the breakdown of superflow in the one- 
band attractive Hubbard model which is a model Hamil- 



tonian for strongly correlated s-wave Fermi superfluids 
on a lattice J^ii^iiiii^ii^ The attractive Hubbard Hamil- 
tonian, 

{ij) 1<T 

i 

describes fermions with spin a hopping, with amplitude 
t, between adjacent sites of a lattice, and interacting via 
an on-site attractive interaction U. The chemical poten- 
tial ^ tunes the fermion density away from half-filling 
at /i = 0. Summation over repeated spin indices is im- 
plicit henceforth, and, unless stated explicitly, we will 
measure energies in units of t. We will study this model 
on the square lattice in two dimensions (2D) and the cu- 
bic lattice in three dimensions (3D). The ground states 
of this model include the uniform superfiuid (SF) and 
a "checkerboard" charge density wave (CDW), which is 
an insulating crystalji^ It is well known that this model 
has an enhanced "pseudospin symmetry" at the point 
/i = 0,>i2, which leads to a degeneracy between the SF 
and the CDW ground states. Away from /x = 0, this 
degeneracy is lifted and the SF ground state is energeti- 
cally more favorable. One expects, however, that so long 
as this pseudospin symmetry is only weakly broken, the 
excitation spectrum of the superfluid will exhibit signa- 
tures of proximity to the CDW ground state. One also 
expects that once superfiow is imposed, there will be a 
critical fiow velocity beyond which the kinetic energy of 
superfiow will overcome the energy difference between the 
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stationary SF and the CDW. The superflow in this case 
is Hmited not by the pairing gap, which may be large, 
but by this energy difference, which becomes very small 
close to half-filling. This situation is similar to the one 
encountered, for example, in cuprate superconductors, 
where the superconducting ground state in the under- 
doped region of the phase diagram is proximate to the 
antiferromagnetic Mott insulator^ Similar physics may 
also be relevant to Cu2;TiSe2 where there is a competition 
between charge order and superconductivityii^ 

This general idea motivates us to study the collective 
modes as well the stability of superflow in the SF phase 
of the attractive Hubbard model. We begin by studying 
the collective mode spectrum in the absence of superflow 
and find that there is indeed a large window of interac- 
tions and density where the collective mode spectrum 
exhibits, in addition to a sound mode, a well-defined 
roton minimum at wavevectors 11 = (tt, tt) (in 2D) or 
n = (tt, tt, tt) (in 3D) arising from proximity to the CDW. 
The existence of such a roton minimum has been pointed 
out in earlier work .^°'^^i^^ We present our results for the 
values of sound velocity and roton gap as functions of 
fermion filling and interaction strength. These results of 
the collective mode spectrum could potentially be veri- 
fied in studies of collective modes in the superfluid phase 
of cold Fermi gases in an optical lattice. We then turn 
to a study of the stability of uniform superflow in the 
Hubbard model. 

We elucidate a variety of superflow-breakdown mech- 
anisms - depairing, Landau instabilities and dynamical 
instabilities. We observe Landau instabilities of the col- 
lective mode at incommensurate wavevectors as has also 
been reported in Ref. [2^ . More importantly, we discover 
dynamical instabilities involving checkerboard or incom- 
mensurate stripe-like density modulations which are dis- 
tinct from previously studied dynamical instabilities of 
Bose superfiuidsi^ii We show that the checkerboard dy- 
namical instability can be viewed as the result of the soft- 
ening of the roton mode and that this instability has a 
simple analog in the strong coupling pseudospin model)^^ 
which we discuss. In contrast, the incommensurate CDW 
instability occurs only at intermediate coupling strengths 
and has no analog in the strong-coupling pseudospin 
model. We explain this instability as a finite momentum 
pairing instability of Bogoliubov quasiparticles, some- 
what analogous to the Halperin-Rice exciton condensate 
instability in indirect band-gap semiconductors.^^ These 
dynamical instabilities could be experimentally probed 
by creating a "running" optical lattice as has been done 
for Bose superfluids.^ 

The existence of the checkerboard dynamical insta- 
bility was discussed earlier by us from the strong cou- 
pling perspectivcj^ The present paper goes consider- 
ably beyond our earlier work. The formalism that we 
use (the generalized random phase approximation or 
GRPA), allows us to address the entire range of inter- 
action strengths from weak to strong coupling. We have 
checked that the GRPA results smoothly match on to 



the earlier analysis of the strong coupling pseudospin 
model and we present these comparisons where appro- 
priate. Within this formalism we have obtained dynami- 
cal phase diagrams of the flowing superfluid over a wider 
range of interactions than in our previous study and un- 
covered the incommensurate CDW instability which was 
not reported in our earlier paper. 

We have tried to make this paper fairly self-contained. 
We begin in Section H with an overview of some key gen- 
eral results on the Hubbard model. Section III focuses on 
the formalism which we use in the remainder of the paper. 
We then turn, in Section IV, to results for the collective 
mode spectrum in the absence of superflow. Section V 
contains the results for the instabilities of the flowing su- 
perfluid and the phase diagram showing the leading insta- 
bilities as a function of interaction and fermion density. 
Section VI contains a discussion of the experimental ob- 
servability of our results in fermionic cold atom systems 
as well as possible implications of these results for vortex 
core physics in strongly correlated superfluids. We end, 
in Section VII, with a summary of our results and possi- 
ble avenues for future research. Appendix A contains the 
Fourier transform convention which we use, Appendix B 
gives details of the GRPA formalism, and the mean field 
theory of the "flowing supersolid" state is discussed in 
Appendix C. 



II. OVERVIEW OF THE MODEL 

Before we turn to a study of collective modes and the 
stability of superflow using the Hubbard model Hamilto- 
nian, we begin by reviewing some well known facts about 
the model which will serve to set the stage for our dis- 
cussion in the remainder of this paper. 



A. Pseudospin Operators and Pseudospin 
Symmetry 

We start by deflning pseudospin operators, first intro- 
duced in the context of the theory of superconductivity 
by Andersoni^ 

rp+ _ t t 

T: = (2) 

where rji = +1 on one sublattice and rji = — 1 on the other 
sublattice of the square or cubic lattice. The physical 
meaning of these operators is clear: creates a fermion 
pair at site i, T~ annihilates a fermion pair at site i, and 
is the pair number operator. It is straightforward to 
check that these operators obey usual spin commutation 
relations. Furthermore, provided that /x = 0, the global 
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pseudospin operators, 

i 

T± = ^7;±, (3) 

i 

all commute with the Hubbard Hamiltonian in Eq. ([1]) so 
that there is a global pseudospin SU(2) symmetry at this 
special point This is in addition to the ordinary global 
spin SU(2) symmetry which is present for any value of ^. 
Pseudospin language is often very convenient and intu- 
itive, in particular when discussing symmetry properties 
and collective modes of paired-fermion superfluids. 



B. Ground States and Degeneracies 

Quantum Monte Carlo simulation a^^i^^'^^ have shown 
that the ground state of the attractive Hubbard model is 
a uniform SF for generic values of U /t and /i/t. The uni- 
form SF has an order parameter {c\-^c\^ ^ Ae^'^ where 
a choice of the phase, ip, of the order parameter corre- 
sponds to a spontaneously broken symmetry. In terms 
of pseudospin operators, this means that (T^) ~ r]iAe^'^ . 
For /i = 0, the pseudospin symmetry discussed above 
leads to the conclusion that all states related to this SF 
ground state by a global pseudospin rotation are also 
valid ground states and are characterized by an order pa- 
rameter N = ?7i(Ti) which can point to any location on 
the Bloch sphere. The uniform SF state has N pointing 
to locations on the equator. The state with N ^ ±Az, 
which points to the north/south pole of the Bloch sphere, 
corresponds to a checkerboard charge density wave state 
(CDW) of pairs. Other locations on the Bloch sphere 
correspond to states with coexisting CDW and SF or- 
ders. 



C. Strong coupling pseudospin Hamiltonian 

In the limit of large U/t, each site on the lattice 
will have either no fermions or a pair of tightly bound 
fermions. Thus, any allowed fermion configuration will 
have Tf = ±1/2 at each site and the kinetic energy will 
lead to tunneling between such "classical" Ising config- 
urations of the pseudospin. The pseudospin dynamics 
for U/t ^ 1 is then described by an effective pseudospin 
Hamiltonian: 

iJeff = J ^ T,.T,~^^T/, (4) 

<ij> i 

where J = At'^/U. For fJ, — 0, this Hamiltonian reduces to 
a pseudospin Heisenberg model which manifestly exhibits 
the pseudospin SU(2) symmetry. 



III. FORMALISM 

The most general problem we would like to tackle, mo- 
tivated by the issues discussed in the introduction, is to 
understand the collective mode spectrum of a flowing su- 
perfluid. We begin by setting up a two-pronged approach 
to attack this problem. At small values of U/t, we build 
upon the mean field theory of the flowing superfluid to 
obtain the collective mode spectrum using GRPA. At 
strong coupling, we study the pseudospin model using 
Holstein-Primakoff spin wave theory. We find that the 
GRPA results smoothly match on to the results from spin 
wave theory of the pseudospin model as we increase U/t, 
showing that the GRPA correctly captures aspects of the 
strong coupling limit. This suggests that the GRPA may 
be a very good approximation to compute the collective 
mode spectrum and address the breakdown of superflow 
over the entire range of couplings. 

A. Generalized Random Phase Approximation for 
the Attractive Hubbard Model 

We begin by constructing the mean field theory of the 
uniform superconducting state of the attractive Hubbard 
model in the presence of nonzero superfiow. The mean 
field Hamiltonian in the presence of a supercurrent is ob- 
tained by forming Cooper pairs with nonzero momentum, 

AoX(e^'^-cl,ct^+e-Q-c,,c,,), (5) 

i 

where we have set U{c^^c^.^) = Aoe**^ ""* and absorbed 
the uniform Hartree shift into the chemical potential. In 
momentum space, the mean field Hamiltonian takes the 
form 

Hmft^YI ^kcLck^-AoX!(4TC-k+Qi + C-k+QiCkT) '(6) 

k k 

where i^k = ~2iek — /i, with ek = X^iLi cos{ki) (c? = 2, 3 
is the dimensionality of the lattice) . 

We can diagonalize Hmft by defining Bogoliubov 
quasiparticles (QPs), 7, via 

f CkT \ _ f Uk(Q) fk(Q) \ f 7kT \ 

I clk+Qx J " I --''(Q) "•^(Q) J \ 7lk+Qx J ■ 

For simplicity of notation, we will refer to the Bogoliubov 
transformation coefficients above as Uk, Wk with the im- 
plicit understanding that they depend on Q. Parametriz- 
ing Uk = cos(0k), I'k = shi{9k), and demanding that the 
transformed Hamiltonian be diagonal leads to the condi- 
tion 

tan(2ek) = TTT^ V 

2 (,4k + 4-k+Qj 
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Defining 



rk = \/^(a + C-k+Q)2 + A2, 



(9) 



wc find that the Bogofiubov transformation coeSicients 
must satisfy the relations 



ul = 




vl = 






Ao 

2rk' 



2rk 



(10) 



In terms of the Bogoliubov QPs, the mean field Hamil- 
tonian finally takes the form 



Hmft = Egs+Y, ^k7L7k<T' 



(11) 



where Egs denotes the ground state energy of -ffMFT and 
denotes the Bogoliubov QP dispersion given by: 



Ek = + 2 (^k - ^-k+q) , 
Egs = ^(a-^^k). 



(12) 



Demanding self-consistency of the mean field theory 
yields the gap and number equations: 

Jj = ^E^(l-"^(^k)-n^(i^-k+Q)), 

k 



f = HME^) + vlil - n^(E_k+Q))] , 



(13) 



where / is the filling, i.e. the average number of fermions 
per site, and A'' is the total number of sites. For given U, 
filling / and fiow momentum Q, these equations can be 
solved to obtain the SF order parameter Aq and the QP 

spectrum. 

Going beyond mean field theory, we need to include 
fluctuations of the density and the superfluid order pa- 
rameter, which we will do within GRPA. We begin by 
considering perturbing external fields that couple to the 
density and order parameter modulation operators as: 



= Hmft -^[hp{i,t)pi 

i 

+ /iA(i,t)Aie*Q--^ +/i^(i,i)Ate-^<^-'-*],(14) 



where 



Ai = CiiCii^. 



(15) 
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FIG. 1: Illustrative example of the collective mode dispersion 
and quasiparticle-pair continuum in 2D, for zero superflow 
(Q ~ 0) with U /t = 3 and / = 0.2 fermions per site, along 
the contour displayed in the inset. 



Going to momentum space, 

i^MFT = Hmft - ^ /i„(K,t)6t (k), (16) 



K 



where a = 1,2,3, summation over a is implicit, and 
6'I'(K) = {/5_K, A_j^_|_Q, A|^_|_q} is the vector of fermion 
bilinear operators corresponding to density and super- 
fluid order parameters at nonzero momenta, with the or- 
der parameters given by: 

Pk = ^1^4 

1 

k 

At = Vrt rt 

— / .'^kt'^-k-l-KI - 
k 

The perturbing flelds correspond to 



^ k 



^-k-l-Ki^kT' 



(17) 



/ii(K,t) = /ip(K,t), 

/l2(K,i) = /lA(K,t), 

/i3(K,i) = /il(-K,i). 



(18) 



We deflne the susceptibility matrix as 

(6„(K,c^)) = x°;3(K,c^)/i^(K,a;). (19) 

The derivation of this bare susceptibility matrix as well 
as the explicit expressions for the elements of this matrix 
are given in Appendix B. 

In order to account for interaction effects at the GRPA 
level, we note that the interaction effectively induces in- 
ternal flelds which renormalize the applied external fleld 
so that the susceptibility within the GRPA can be ex- 
pressed as: 



Xai3 



(20) 



Here, D = diag{2, 1,1} is a diagonal matrix that en- 
codes the decoupling of the interaction Hamiltonian. The 



derivation of this expression in explained in Appendix B. 
This GRPA susceptibility will diverge when the determi- 
nant Det(l — Uy^D) becomes zero (or equivalently, one 
of the eigenvalues of this matrix vanishes). We identify 
the locus of real frequencies, uj = ^^'(K), where this hap- 
pens, as the dispersion of a sharp (undamped) collective 
mode. 

Fig. [1] provides an illustrative example of the collective 
mode spectrum obtained using the GRPA in 2D with 
U /t = 3 and a filling / — 0.2 fermions per site in the 
absence of superflow {Q — 0). We find a linearly dispers- 
ing superfluid "phonon" mode at small momenta and low 
energy. The slope of this linear dispersion is the sound 
velocity. The collective mode disperses as a function of 
momentum over the Brillouin zone and exhibits an ex- 
tremum at the edges of the Brillouin zone. In the example 
here, it is a maximum at K = (tt, tt), but with increas- 
ing interaction strength and at fillings closer to / = 1, 
the spectrum exhibits a minimum at K = (tt, tt) arising 
from strong short distance density correlations. At high 
energies, there is an onset of a two-quasiparticle contin- 
uum where the collective mode can decay by creating two 
Bogoliubov quasiparticles with opposite spins in a man- 
ner which conserves energy and momentum. Once the 
collective mode energy goes above the lower edge of the 
two-particle continuum of Bogoliubov QP excitations, it 
will cease to be a sharp excitation and acquire a finite 
lifetime. The collective mode energy as well as the pair 
continuum change in the presence of superflow as dis- 
cussed subsequently. 



B. Strong Coupling Limit: Spin Wave Analysis of 
Pseudospin Model 




FIG. 2: (Color online) Collective mode energy at zero super- 
flow in 2D at strong coupling, U/t = 15. The dispersion is 
shown along the indicated contour in the Brillouin zone for 
different fillings: (a) / — 0.8 fermions per site and (b) / = 1.0 
per site. The GRPA result (solid line) is in good agreement 
with the Holstein-Primakoff spin wave result (dashed line, 
HP) for the strong coupling pseudospin model. The roton 
minimum has a small gap at / = 0.8 but becomes a gap- 
less mode at / = 1.0 due to the pseudospin SU(2) symmetry 
discussed in the text. 



In the pseudospin model, a state with nonzero super- 
current is obtained by imposing a phase twist on the 
non-flowing ground state |0) as 



IQ) = exp 



|0). 



(21) 



Equivalently, we can make a unitary transformation to 
work with the Hamiltonian 



i?eff(Q) 



(I^fTf + TfTf)cos(Q-r, 



(j:-T/-I^«T;)sin(Q.r,,)] 



(22) 



where = r,; 



Tj. This amounts to transforming to a 
reference frame where the superfluid is at rest. 

The classical ground state |0)c, which is the ground 
state obtained under the assumption that the pseu- 
dospins are classical vectors, is equivalent to the mean- 
fleld ground state of the interacting fermion model and 



can be parametrized by specifying the classical vector 
at each lattice site. This is given by 



= S{T]i sin6l,0,cos6'). 



(23) 



where the pseudospin magnitude S — 1/2. The angle 9 
is related to the filling, /, defined earlier, as 



/ — 1 — cos ( 



(24) 



Working at fixed filling, the chemical potential is given 
by 



^l = 2JScos{9){eo + eg), 



(25) 



where eg = J2i=i cos{Qi) as before. 

We compute the excitations about the ground state 
|0)c for the Hamiltonian i7cff(Q) using a Holstein- 
Primakoff (HP) approach^^ by rotating to new pseu- 
dospin operators T given by 



f: = Tlcos{9)+i^,T:sin{9), 
f[ = -T:sm{9)+v,T,^cos{9), 



(26) 



and expressing T in terms of HP bosons. Expanding the 
Hamiltonian to 0(5*), we find 

H = E, + 5E^ + Y^u;^{Cl)hy)^, (27) 



K 



where 



^c==-iVJS'^[cos^ Oea + (f + cos^ 6')eQ], 
JS_ 
T 



JSm-^, . 2n (l+cos^fl), 



K 



represent, respectively, the classical ground state energy 
and the leading quantum correction to the ground state 
energy. The spin- wave dispersion clIk(Q) is given by 



u;k(Q) = 2J^(/3k(Q) + V«k(Q)-7^(Q)), (29) 



with 



rr^\ silled (l+cos2 6') /eK+Q + EK-Q 
aK(Q) =eQ + ^:^eK ^ 



^k(Q) = -cos6l(eK-Q - eK+q) , 

7k(Q) = Jsin2 0(2eK + eK+Q + eK-Q). (30) 

An illustration of the collective mode dispersions ob- 
tained using this approach is shown in Fig. [2] at a strong 
coupling value of the interaction U/t — 15 for two dif- 
ferent fillings, / = 0.8, 1.0 fermions per site. We find 
that these dispersions are in good agreement with those 
obtained using the GRPA which serves to show that the 
GRPA also correctly captures aspects of the strong cou- 
pling limit. As seen from the dispersion, there is a lin- 
early dispersing superfluid "phonon" mode at small K. 
In addition, for / = 0.8, we find a low energy "roton min- 
imum" a,t K = (tt, tt) which arises from strong local pair 
density correlations in the superfluid ground state. For 
/ — 1.0, the gapless roton mode at K — {tt, tt) is a man- 
ifestation of the pseudospin SU(2) degeneracy, discussed 
earlier, between the crystalline CDW ground state and 
the uniform superfluid ground state. 
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FIG. 3: (Color online) The sound mode velocity, Vs, as a 
function of the fermion filling / in (a) 2D and (b) 3D for 
U/t = 3 (circles) and U/t = 15 (diamonds). Solid line is 
a guide to the eye. The dashed lines are the weak coupling 
result, Vs = (^;F/\/d)[l-^7iV(0)]^/^ from Ref.^ for U/t = 3, 
with A''(0) being the non-interacting density of states (per 
spin) at the Fermi level. The dotted line indicates the Holstein 
PrimakofF spin- wave result for U/t — 15. The inset to (b) 
shows the expected t/U scaling of Vs/t for U/t ^ 1 in 3D. 



rotons only exist at this commensurate wave vector on 
the lattice and do not form an incommensurate ring of 
wave vectors. Further, the roton mode becomes gapless 
at /X = and this gaplessness arises from a pseudospin 
SU(2) symmetry which leads to a degeneracy between the 
superfluid and CDW ground states as discussed earlier. 



IV. COLLECTIVE MODES IN THE 
STATIONARY SUPERFLUID 

We have seen that the collective mode spectrum ob- 
tained using the GRPA or the strong coupling analysis 
has a superfluid phonon mode at small momenta and a 
roton minimum at K = (tt, tt) in 2D or at K = (tt, tt, tt) in 
3D. The phonon mode is a reflection of the fact that the 
superfluid ground state breaks a global U(l) symmetry 
- it is thus the long wavelength Goldstone mode asso- 
ciated with this symmetry breaking. The roton mode 
arises from strong local density correlations, analogous 
to the roton mode in ^He. However, unlike in ''He, the 



A. Sound Velocity 

Fig. [3] shows the sound mode velocity Vs, extracted 
using the GRPA, over a range of densities and interac- 
tion strengths in 2D and 3D. In the limit of low filling 
and weak interaction, our calculations of Vs agree with 
the results of Belkhir and Randeria,^- although strong 
finite size effects and a small pairing gap limit our nu- 
merical exploration of the very weak coupling regime 
U/t <C 1. At strong coupling, U/t^ 1, our results from 
the GRPA are in good agreement with the spin wave re- 
sult = {At'^/U)Vd^2f- p. The linear scaling of Ws/i 
with t/U , expected for C//i ^ 1, is illustrated in the inset 
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FIG. 4: (Color online) The energy of the collective mode at 
(a) (tTjTt) in 2D and (b) (7r,7r,7r) in 3D for different inter- 
action strengths. The dashed (solid) lines indicate that the 
mode energy corresponds to a local maximum (minimum) of 
the dispersion. The inset shows the collective mode energy in 
3D at (tt, tt, tt) at a filling of / = 0.8 fermions per site, as a 
function of t/U . Inset shows a comparison of the GRPA re- 
sult (crosses) with the strong coupling spin-wave theory result 
(dashed line). 



of Fig.[3);b). 



B. Roton Gap 

Fig. [4] shows the energy of the collective mode dX K — 
(tt, tt) in 2D and K = (7r,7r,7r) in 3D. For low fillings 
and weak to intermediate coupling strengths, the mode 
energy shows a maximum at this momentum. For strong 
enough interactions, or for fillings near f — 1.0, however, 
the mode energy exhibits a local minimum and can be 
justifiably identified as a roton mode. The roton gap 
clearly scales as t/U for U/t ^ 1 as seen from the inset 
of Fig.|4|D. Further, the roton energy goes to zero linearly 
as / — > 1. This is due to the chemical potential /i, which 
tunes the filling / away from half-filling, and explicitly 
breaks the pseudospin SU(2) symmetry present at / = 1. 



V. SUPERFLOW INSTABILITIES 

Superflow instabilities are easily understood in the con- 
text of bosonic superfluids possessing Galilean invariance. 
Superflow can then be simply thought of as a Galilean 
boost performed on a stationary SF. The excitations in 
the flowing frame are the same as that of the stationary 
SF, except their energies undergo a Doppler shift. At 
a certain critical flow velocity, the Doppler shift renders 
the excitations gapless at some momentum. These gap- 
less excitations are then populated, flow energy is trans- 
ferred to these bosonic excitations, and superfluidity is 
lost. The critical velocity is given by Wcrit = iiiiii(wg/'j'), 
where LOq is the energy of an excitation of the stationary 
SF at momentum q. 

Our system, a paired fermionic superfluid on a lattice, 
is much richer. The excitations are fermionic Bogoliubov 
quasiparticles and a bosonic collective mode, correspond- 
ing to coupled quantum-mechanical fluctuations of the 
superfluid order parameter and the density. The effect 
of flow is now more than simply introducing a Doppler 
shift in excitation energies since the system explicitly 
breaks Galilean invariance. The SF order parameter and 
the quasiparticle dispersion of the mean field theory are 
strongly renormalized. In addition, the collective mode 
dispersion also changes strongly. This is easily seen in 
the strong coupling pseudospin model, where the disper- 
sion of spin wave collective modes is given by Eg. ip^)) . 
The term /3k (Q) in this equation is the Doppler shift 
and clearly depends on the flow momentum. However, 
the terms aK(Q) and 7k (Q) are strong functions of the 
flow momentum too. Further, on the lattice as opposed 
to the continuum, the Doppler shift vanishes at momen- 
tum points corresponding to the Brillouin zone edges. 
The only effect of superflow in this case is to cause a 
strong dispersion renormalization via its effect, at strong 
coupling, on aK(Q),7K(Q)- 

Due to these effects of imposed superflow, our sys- 
tem undergoes three broad kinds of instabilities, which 
we call "depairing" , "Landau" and "dynamical" . Be- 
low, we describe each of these and map out a stability 
"phase diagram" as a function of dimensionality, inter- 
action strength and density, indicating the critical flow 
momentum beyond which the uniform flowing superfluid 
is unstable. Here we will restrict our attention to super- 
flow momenta Q = Q^x. 



A. Depairing Instability 

The system undergoes a depairing instability when the 
self-consistently calculated SF order parameter. A, van- 
ishes in the mean field theory of the flowing SF. At the 
critical fiow momentum, the energy cost of flow outweighs 
the condensation energy gain, and the system goes into 
an unpaired normal state. 

This instability is close to, but not identical with, the 
quasiparticles becoming gapless due to the flow-induced 
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Doppler shift. In 2D, we find that these two phenomena 
occur at the same value of the flow momentum. In 3D 
however, superfluidity persists beyond the point where 
quasiparticles become gapless. There is a small window 
of gapless superfluidity, where negative energy Bogoli- 
ubov quasiparticle states are occuppied. This is analo- 
gous to what has been claimed to occur, for instance, in 
superfluid ■^He in the presence of superflow.^^ In this case, 
although there are negative energy quasiparticle excita- 
tions, the system cannot arbitrarily lower its energy by 
occupying such states since the Bogoliubov excitations 
arc fcrmionic. The gapless superfluid thus survives as a 
stable intermediate phase in 3Di^ 



B. Landau Instability 

A Landau instability occurs when the collective mode 
energy hits zero and becomes negative, as shown in Fig.O 
In the strong coupling limit, with the collective mode 
dispersion given in Eq.(29), this case is described by 

«k(Q) > 7k(Q) and /3k (Q) < - V"k(Q) + 7k(Q)- 
In the GRPA approach, we find that the renormalized 
susceptibility, x'^^^^: diverges for two real negative fre- 
quencies. These frequencies match smoothly onto the 
results of the spin wave frequencies at strong coupling. 

As discussed in Ref.|23[, Landau instability is not an 
instability of linearized dynamics of small fluctuations 
around the uniformly flowing state, since some form of 
mode coupling is necessary to transfer the energy of the 
superflow into these negative-energy modes. Its full the- 
oretical description is thus rather complicated. Here we 
restrict ourselves to only finding the critical flow momen- 
tum for this instability. 

We find that the Landau instability occurs either at 
small momenta, corresponding to the dispersion of the 
sound mode going below zero, or at some finite incom- 
mensurate momentum. In the case of low filling, unless 
preempted by depairing, we see a Landau instability of 
the sound mode. For moderate values of U and the fill- 
ing, we see Landau instabilities at large incommensurate 
momenta (as in Fig. [5]). As the filling is reduced, this 
incommensurate wavevector moves towards the Brillouin 
zone centre, so that in the low density limit, it is the long 
wavelength sound mode that becomes unstable. Simi- 
lar incommensurate instabilities have also been found in 
Ref.fH. 



C. Dynamical Instability 

The underlying lattice potential allows for a new pos- 
sibility for the breakdown of superfiow, namely through 
dynamical instabilities.^ This corresponds to the eigen- 
frequency of the collective mode dispersion in our system 
becoming zero, and subsequently becoming complex. 

The concept of dynamical instability is particularly 
simple to understand for weakly interacting bosons. In 
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FIG. 5: Landau instability - Collective mode spectrum (from 
GRPA) in 2D with U jt = 7 at a filling of / = 0.4 fermions per 
site. Collective mode dispersion is shown at (a) zero flow, (b) 
just below, and (c) just above, the critical flow momentum 
for the Landau instability. The -|-(— ) sign indicates the mode 
is at wavevector -|-K(— K) which has an s-component along 
(opposite to) the flow direction so that it is Doppler shifted 
down (up) in energy. The collective mode frequency becomes 
negative at an incommensurate wavevector. 



this case, the instability coincides with the point where 
the effective mass of the bosons changes sign as a func- 
tion of the superflow momentum,— leading to runaway 
growth of phase and density fluctuations which eventu- 
ally destroy superfluidity. 

For strongly interacting bosons the situation is more 
interesting. It has been shown^ that with increasing 
interaction strength at a commensurate flUing (integer 
number of bosons per site), the dynamical instability oc- 
curs at a smaller and smaller superflow momentum. The 
critical flow momentum eventually tends to zero at the 
equilibrium superfluid to Mott insulator transition, scal- 
ing as the inverse of the divergent correlation length as- 
sociated with this Mott transition. 

Here, in the fermionic analog of the problem of the 
critical superflow in a strongly correlated superfluid, we 
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FIG. 6: Dynamical Instability (Commensurate)- Collective 
mode spectrum (from GRPA) for U/t = 5, filling / = 0.8 
fermions per site on a 2D square lattice. The +(— ) sign in- 
dicates the mode is at wavevector -|-K(— K) which has an 
a;-component along (opposite to) the flow direction. As the 
flow momentum Q is increased, the collective mode frequency 
at K = (7r,7r) decreases until it hits zero and becomes com- 
plex. This gives rise to a dynamical instability associated with 
the "checkerboard" CDW order. The part of the dispersion, 
around (ti,h), which corresponds to unstable modes is not 
shown. 



find, in addition, an entirely different kind of dynamical 
instability, associated with emergent density-wave order 
for a large range of fillings, wherein the mode energies be- 
come complex at nonzero wavevectors. In the strong cou- 
pling limit, this happens when |q!k(Q)| becomes smaller 
than |7k(Q)|- In the GRPA approach, at a dynamical 
instability, we find no real frequency at which any eigen- 
value of the inverse GRPA susceptibility vanishes. The 
GRPA results for the critical flow momentum Q, and the 
wavevector K at which the dynamical instability occurs, 
match smoothly onto the spin wave results at strong cou- 
pling. In this limit it is easy to check analytically that 
the dynamical instability associated with this emergent 
charge order is not preempted by a Landau instability. 



The appearance of complex collective mode frequencies 
indicates an exponential growth of density fluctuations 
around the uniformly flowing state, at a wavevector K. 

We find two kinds of such dynamical instabilities — 
associated with commensurate (Fig. 15]) or incommen- 
surate (Fig. [7]) charge order. We see the former for a 
large range of fillings in the vicinity of / = 1, for all 
interaction strengths. The wavevector at which this in- 
stability happens is at the Brillouin zone corner - (tt, tt) 
in 2D and (tt, tt, tt) in 3D. This corresponds to an instabil- 
ity towards a "checkerboard" CDW state, with a density 
modulation of opposite sign on the two sublattices of the 
square or cubic lattice. The incommensurate dynamical 
instability, on the other hand, occurs at wavevectors cor- 
responding to various incommensurate ordering patterns 
and it arises as follows. In the presence of superfow, the 
dispersion of the Bogoliubov quasiparticles is renormal- 
ized and Doppler shifted, giving rise to multiple minima 
separated by a nonzero wavevector. This leads to a peak 
in the bare mean field susceptibility, at this momen- 
tum. The interaction renormalizes this peak into a diver- 
gence and leads to a finite momentum pairing instability 
of Bogoliubov quasiparticles. This dynamical incommen- 
surate instability is thus a non-trivial, interaction-driven 
phenomenon, which only occurs at intermediate coupling 
strength. This phenomenon is somewhat analogous to 
the formation of cxcitonic condensates in indirect band- 
gap semiconductors)^ 



D. Stability phase diagrams 

Taking these instabilities into account, we map out su- 
perfiow stability phase diagrams in 2D (Fig. [8]) and 3D 
(Fig. [9]). These plots show the first instability that is 
encountered as imposed fiow is increased, for different 
values of filling. Values of U for the plots have been cho- 
sen so as to illustrate all possibilities. We see that the 
commensurate dynamical "checkerboard" CDW instabil- 
ity comes into play around / = 1/2 for all values of the 
interaction strength and is the dominant instability all 
the way to / = 1, where the critical flow momentum 
vanishes, reflecting the degeneracy between the SF and 
CDW states. 

In the low density limit, the system is similar to a 
continuum Fermi gas. In 3D, the density of states at 
the Fermi level vanishes. At low interaction strength, 
this leads to the pairing gap A being exponentially sup- 
pressed. The sound velocity, on the other hand, is pro- 
portional to the Fermi velocity, so that Vs f^^^- This 
leads to a rather sharp drop of the sound velocity as 
/ ^ but the gap drops to zero much faster. There- 
fore, in the 3D case, at weak interaction and low filling, 
a small imposed flow will drive A to zero before the flow 
velocity exceeds the sound velocity, leading to a depairing 
instability as can be be seen in Fig. [51 

In contrast, in the low density continuum limit in 2D, 
the density of states goes to a constant, which means the 
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FIG. 7: Dynamical Instability (Incommensurate) - Collective 
mode spectrum (from GRPA) for U/t — 3, filling / = 0.6 
fermions per site on a 2D square lattice. The +(— ) sign in- 
dicates the mode is at wavevector -|-K(— K) which has an x- 
component along (opposite to) the flow direction. As the flow 
momentum is increased beyond 0.27r, the collective mode 
frequency becomes complex at an incommensurate wavevec- 
tor. The complex frequencies in the unstable region are not 
shown in the figure. 
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FIG. 8: (Color online) Stability phase diagram for a 2D square 
lattice case with (a) U/t = 3 and (b) U/t = 5. For every fill- 
ing, the plot shows the first instability that is encountered as 
the fiow is increased. The solid (blue) line in the low den- 
sity limit in (a) indicates the region where we expect to see 
a Landau instability, but finite size effects prevent us from 
accessing the area. The other transitions in the figure cor- 
respond to Depairing (circles, dotted line), Incommensurate 
Dynamical Instability (diamonds, solid line), Commensurate 
Dynamical Instability (squares, dashed line), Landau Insta- 
bility (triangles, dash-dotted line). At very large interaction, 
the diagram looks similar to U/t — 5, except that the incom- 
mensurate dynamical instabilities disappear. 



pairing gap stays finite as / — > 0. With imposed flow, 
the first instability that one encounters is then the Lan- 
dau instability, which will happen when the flow velocity 
exceeds the sound velocity which scales as Vs /^^^. We 
have not been able to numerically uncover the Landau 
instability in this regime due to severe finite size effects. 

At intermediate densities, and at intermediate values 
of the interaction strength, the leading instability is the 
incommensurate dynamical instability, corresponding to 
the emergence of an incommensurate CDW, character- 
ized by a wavevector which depends on the density and 
the interaction strength. This instability does not hap- 
pen in either the strong-coupling limit, where the Ander- 
son pseudospin description is appropriate^ nor in the 
weak-coupling BCS superfluid limit, where the instabil- 
ity at intermediate densities is due to depairing. This dy- 



namical instability is a nontrivial intermediate-coupling 
phenomenon and is one of the most interesting results 
reported in this work. 

An interesting question that arises in relation to the 
emergent CDW-driven dynamical instabilities, is the fate 
of the system beyond the critical fiow momentum. We 
show, in Appendix C, that beyond the commensurate 
CDW dynamical instability threshold the system behaves 
as a "flowing supersolid" at the mean field level. How- 
ever, we argue that fluctuations beyond mean field theory 
are expected to destroy this flowing supersolid. An ex- 
plicit example of this, in the strong coupling limit, has 
been discussed by us earlier 3. We expect the same to be 
true in the case of the incommensurate CDW dynamical 
instability. 
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FIG. 9: (Color online) Stability phase diagram for a 3D cubic 
lattice with (a) U/t = 5 and {h)U/t — 15. For each filling, 
the instability first encountered as the flow is increased is 
shown. The transitions in the figure correspond to Depairing 
(circles, dotted line). Incommensurate Dynamical Instability 
(diamonds, solid line). Commensurate Dynamical Instability 
(squares, dashed line). Landau Instability (triangles, dash- 
dotted line). 



VI. EXPERIMENTAL IMPLICATIONS 



moving optical lattice. In this frame, the bosons must 
condense in a state of non-zero lattice momentum. If the 
optical lattice is moved sufficiently slowly, we expect the 
initial zero-momentum SF ground state of the bosons to 
adiabatically acquire a phase gradient, thus increasing 
the superflow momentum until dissipation sets in and 
destroys the condensate upon crossing a Landau or dy- 
namical instability. 

The instabilities of superflow for fermions loaded in an 
optical lattice may similarly be tested in cold atom exper- 
iments. This has been clearly demonstrated in a pioneer- 
ing experiment of Ref. [TT| which studied the superfluid 
phase of fermions in a harmonic trap and showed that 
the critical velocity smoothly interpolates between the 
pair breaking velocity in the BCS regime to the Landau 
critical velocity in the BEC regime.-" In a deep optical 
lattice, provided that the superflow does not lead to exci- 
tations into the higher bands, the superfluid phase should 
be reasonably well described by the Hubbard model in 
the lowest band. In this case, one may be able to probe 
the various dynamical phase diagrams proposed in this 
paper. Even before reaching the instability, one may be 
able to use Bragg spectroscopj*2i to probe the collective 
density fluctuation spectrum, and observe the roton soft- 
ening with increasing lattice velocity which would enable 
one to test the mechanism by which the superflow breaks 
down. 

Our results are also of relevance for the vortex-core 
physics in strongly-correlated lattice superfluids. Con- 
sidering the vortex core as simply the region where the 
current has exceeded the value beyond which the uniform 
superfluid phase is no longer stable, we would expect the 
dynamical instabilities discussed in this paper to lead to 
commensurate or incommensurate density orders in the 
vortex core. 



Assuming that future experiments will be able to pro- 
duce attractive cold Fermi gases in the lowest band of an 
optical lattice, it is reasonable to expect that the ground 
state and low energy excitations of such a system will be 
well described by the attractive Hubbard modeli^ii^ Our 
results on the collective mode excitations of the Hubbard 
model, specifically the roton gap and the sound velocity 
as a function of filling and interaction strength, could 
then be verified experimentally. 

Turning to superflow instabilities, the theoretically 
predicted Landau and dynamical instabilities of Bose 
superfluids^ii^ have been experimentally observed us- 
ing cold bosonic atoms in an optical lattice.^ As shown 
there, a flow in a cold atom system can be induced by 
frequency-detuning one pair of the counterpropagating 
laser beams that form the lattice. This creates a "run- 
ning" optical lattice potential in the corresponding basis 
direction, i.e. the lattice is effectively moving with re- 
spect to the stationary potential of the trap, holding the 
atomic cloud. The effect of this on the bosons is easy 
to understand if one goes to the reference frame of the 



VII. CONCLUSIONS 

In conclusion, we have studied the SF ground state 
of the negative-U Hubbard model on a square and cu- 
bic lattice. Using GRPA, we have studied the collective 
mode spectrum at all values of the interaction strength 
and density. We have shown that strong density correla- 
tions at increasing filling lead to the appearance of roton 
features in the collective mode spectrum, which signal 
developing competition between uniform SF and density- 
ordered CDW states. When superflow is imposed on the 
system, depairing and Landau instabilities arise. In ad- 
dition, this competition leads to dynamical instabilities 
corresponding to roton mode softening. This dynamical 
instability occurs at a commensurate wavevector corre- 
sponding to the corner of the Brillouin zone. This in- 
stability is easily described in the strong coupling pseu- 
dospin language and is a consequence of the degener- 
acy between SF and "checkerboard" CDW ground states 
at / = 1. Another very interesting dynamical instabil- 
ity occurs at an incommensurate wavevector and has no 
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strong coupling analog. It occurs over a range of densities 
near f = 1/2 and at intermediate values of the interac- 
tion strength. Our results may be verified in experiments 
with ultracold fermions loaded on an optical lattice. In- 
teresting avenues to explore in the future would include 
extending this work to other lattice geometries^ and to 
study superflow instabilities in multiband superfluidsj^ 
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APPENDIX A: FOURIER TRANSFORM 
CONVENTIONS 



where the perturbation is assumed to be Hermitian. Us- 
ing time-dependent perturbation theory to leading order 
in the perturbing fields, one finds that the change in the 
expectation value of the operator Oq,(K) is given by: 



^(Oa)(K,i) 



where 



+ 00 



dt' xlp{^.t-t')hp{V.,t'), (B2) 



X° .(K, t-t') =^0(i-i')( [O JK, i), 6t (K, t')l )„. (B3) 



Here [.,.] denotes the commutator and (.)o implies that 
the expectation value is taken in the ground state of Hq. 
Upon doing a spectral decomposition, one obtains: 



We have Fourier transformed fermion operators as 

^J^Ck.e^''-. (Al) 



This defines the Fourier transform relations of our mod- 
ulation operators as 



K 



(A2) 

(A3) 
(A4) 



where the momentum space operators pK and Ak are 
given in Ea. p7|) . 

Also, the real space fields that couple to our modula- 
tion operators are Fourier transformed as 



(A5) 



These Fourier transforms have been used in going from 
Eq.lUll) to Eq.®. 



Here {0)mn = {m\0\n), and |7i),|m) denote the eigen- 
states of Ho (with n — corresponding to the ground 
state). In the denominator, Eno = En—Eo where £"„ is 
the energy of state \n). 

We are interested in the case where Hq — Hmft, with 
the operators 0}y.(K) given by 



Ol(K) 

oUk) 
oUk) 



(B5) 



and the perturbing fields corresponding to 
hi{K,t) = hp{K,t), 

h2(K,t) - hA{K,t), 

/i3(K,t) = /il(-K,t). 



(B6) 



APPENDIX B: DETAILS OF THE 
GENERALIZED RANDOM PHASE 
APPROXIMATION 

Let us assume that we have a Hamiltonian Hq which is 
modified by a set of weak time-dependent perturbations 
so that 

H(t) = iJo-^^/ia(K,i)Ot(K), (Bl) 

K 



At zero temperature, it is easily checked that the quasi- 
particle energies are all positive in the dynamically stable 
region of interest. In this case, the only terms which have 
nonzero matrix elements in the expression for x^^(K, uj) 
are those where the operators create two Bogoliubov 
quasiparticles (QPs) when acting on |0) and where they 
destroy two Bogoliubov QPs when acting on |n). It is 
therefore convenient to resolve each perturbation opera- 
tor into two parts - one that creates two QPs and one 
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P-K = 



that annihilates two QPs - as follows: 

- ^X!("k+KWk + Ukt^k+K)7k+KT'^-l. 
k 

- ^Xl("k-Kfk + MkUk-K)7-k+Qi7k-KT' 
k 

A1k+Q = -X!^'k«k+K7k+KT7-k+Qi ' 
^K+Q = I]"kWk+K7k+KT7-k+Qi ' 



^-K+Q = X!"kWk-K7-k+Qi7k-KT' 
k 



-Xl^kt^k-K7_k+Qi7k-KT- 
k 



where the superscript 'c' denotes creation of 2 quasipar- 
ticles, and 'a' denotes annihilation. 

Since is a symmetric matrix, it suffices to use the 
above to compute the following distinct elements: 



Al,l 



Al,2 



Al,3 



A2,2 



X3,3 



-E 

AN ^ 



-E 

2iV ^ 



-E 



-E 



-E 



-E 



(Uk-K^'k + t'k-KMk)^ 
+ -Bk-K + £^-k+Q 

_ (ttkfk+K + «kMk+K)'^ 
W — i?k+K — -E_k+Q 
(Uk-K^^k + t'k-KMk)MkMk-K 
+ -Bk-K + £^-k+Q 

(Mfeffe+if-K!fcMfe+if)f/ct^fc+if 



— -Ek+K — -£^-k+Q 
(Mk-Kt^k + t'k-KMk)t^kt^k-K 
+ £^k-K + -B_k+Q 

(MkWk+K+VkWk+K)MkWk+K 



W— ii'k+K— ^'-k+Q 
«k"k-K 



k+Q 



2, ,2 



k'^k+K 



— -Ek+K — -B-k+Q 
Mk^k^k-Kt'k-K 

^ + -Ek-K + £^-k+Q 

^ ttk^kMk+K'^k+K 

— -Ek+K — -E-k+Q 

^k^k-K 
+ -Ek-K + -E_k+Q 

"k^k+K 

W — i?k+K — £^-k+Q 



(B8) 



In order to include interaction effects within the 
GRPA, we note that the interaction can be decomposed 



as follows: 

— Ucl-^cl^CiiCi]' 



(4t 4i > c»i t+(c,x c,| ) (B9) 



These expectation values, generated by interactions in 
the presence of external fields, act as "internal fields" 
which renormalize the applied field. It is easy to show 
that this effect can be taken into account by simply set- 
ting 



This leads to: 



hi{K,Lo) + 2U(di(K,u)), 

h2{K,Lu)+U{d2{K,Lu)), 

/i3(K,c^) + f/(03(K,c^)). 



(BIO) 



(Bll) 



where D = Diag{2,l,l} is a diagonal matrix, and we 
have suppressed (K,aj) labels for notational simplicity. 
Solving the above equation gives: 

<5(0„) = [(1 - Ux°DrW0 hp ^ (B12) 
where x'^p^^ is the renormalized GRPA susceptibility. 



APPENDIX C: MEAN FIELD THEORY OF THE 
FLOWING SUPERSOLID 

For all values of the interaction strength and for 
fermion densities near / = 1, we clearly see a dy- 
namical instability at the commensurate "checkerboard" 
wavevector - (7r,7r) in 2D and (7r,7r, tt) in 3D. At this in- 
stability, we expect "checkerboard" density order to arise 
beyond critical flow, leading to a state with coexisting SF 
and density order, a "flowing supersolid" i^^*^ We study 
the mean field theory of this state in order to examine its 
stability. In the presence of a nonzero SF fiow, imposed 
as a uniform phase gradient on A, we self-consistently 
calculate the ground state allowing for "checkerboard" 
density modulations. 

The mean field order parameters are: 



A = 

P = 
A = 



■^E^^-k+Qi^kr) 



2^ E^('"k+nCT'"ker) 
k 



u 

N 



y^(c-k+n+QiCkT) 



(CI) 



where 11 = (7r,7r) in 2D or (tt, 7r,7r) in 3D. Since there 
are no spin-selective Zeeman terms, we are justified in 
forcing the density modulation to be equal for both spin 
species. Due to global phase rotation U(l) symmetry, we 
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can choose A to be real but we allow A to be complex. 

being the expectation value of the staggered density, 
is real. 

Upto an unimportant constant, the Hamiltonian can 
then be written as: 



H 



k 



(C2) 



where the primed summation in the Hamiltonian indi- 
cates that if k is included, then k + 11 is to be excluded. 
The other notation is as follows: 



*k = 



k+Qi 

Ck+nt 
-k-n+Qi 



and 



H = 



( a 

-A 



V 



-p 



-A -p 

-e-k+Q -A* 

-A a+n 

p -A 



-A 

P 
-A 

C-k-n+Q / 



(C3) 



(C4) 



For given Q and the density, we numerically diago- 
nalize this matrix and solve the self-consistency equations 
for A, p, A and the filling /. We then evaluate the su- 
perfiuid and density order parameters in the converged 
solution. In addition, we also evaluate the uniform cur- 
rent, 



{J) = -2t{Y^clc, 



,Vkek) 



(C5) 



where — 2tek is the non-interacting fermion dispersion. 

A "flowing supersolid" phase is indicated by simultane- 
ous non-zero values for A and p. For the case of checker- 
board order we do see a "flowing supersolid" phase be- 
yond a critical flow momentum. However, the onset of 
density order coincides with a maximum in the expecta- 
tion value of the current as a function of flow momen- 
tumfFig. [TU)) . This indicates that the system is dynami- 
cally unstable to the exponential growth phase and den- 
sity fluctuations. The argument below demonstrates why 
this is the case in a simpler context. 

Let us consider a superfluid system in ID for simplicity. 
Denoting the mean field density and phase of the SF 
order parameter by no and 0o respectively, we consider 
fluctuations Sn and 64>. Now, the current is some function 
of the gradient of the phase: 



(C6) 



The equations governing the dynamics of the fluctuations 
are the Josephson relation and the continuity equation. 
The former gives: 



3 
2.5 
2 



1.5 
1 

0.5 



A 




P.; 






A _ 









0.1 0.2 0.3 0.4 

Flow momentum Q (in units of Jt) 



0.5 



FIG. 10: (Color online) Mean field theory results for the 
"flowing supersolid" state showing the various observables as 
functions of the flow momentum in 2D for U/t=7 and with a 
filling of 0.8 fermions per site. Supersolid order onsets around 
Qx « 0.2n, as both A and p have simultaneous non-zero 
expectation values. This coincides with a maximum in the 
current as a function of the flow momentum, indicating a dy- 
namical instability. 



where a — dp / dn > is closely related to the compress- 
ibility. The continuity equation is: 



dn 



Substituting n — fiQ -\- 
Q=^,we obtain 

dSn 
~dF 



dj 
dx 

5n and 



dx 



Q 



^, where 



dJ d^ 



dQ dx^ 



(C8) 



Combining this with the Josephson relation, we finally 
obtain: 



1 d^ 



dt^ 



dJ d'^6(l) 
dQ dx'^ 



(C9) 



When dJ /dQ, becomes negative, the wavelike solutions 
of this equation develop complex frequencies. In this 
case, density and phase fluctuations will grow exponen- 
tially in time making the system dynamically unstable 
when the current goes through a maximum as a function 
of the flow momentum Q. 
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